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Abstract. This paper deals with the computation of polytopic invariant sets for polynomial dynam- 
ical systems. An invariant set of a dynamical system is a subset of the state space such that if the 
state of the system belongs to the set at a given instant, it will remain in the set forever in the future. 
Polytopic invariants for polynomial systems can be verified by solving a set of optimization problems 
involving multivariate polynomials on bounded polytopes. Using the blossoming principle together 
with properties of multi-afiine functions on rectangles and Lagrangian duality, we show that certified 
lower bounds of the optimal values of such optimization problems can be computed effectively using 
linear programs. This allows us to propose a method based on linear programming for verifying 
polytopic invariant sets of polynomial dynamical systems. Additionally, using sensitivity analysis of 
linear programs, one can iteratively compute a polytopic invariant set. Finally, we show using a set 
of examples borrowed from biological applications, that our approach is effective in practice. 



An invariant set of a dynamical system is a subset of the state space such that if the state of the 
system belongs to the set at a given instant, it will remain in the set forever in the future. Invariant 
sets are fundamental notions in dynamical systems theory where they can serve to prove the existence 
of attractors (e.g. in the Poincare Bendixon theorem |Sas99j ): they have played an important role 
in control theory for the analysis of performance, robustness or practical stability [Bla99]. They 
are also of great interest for reachability analysis of continuous and hybrid systems, especially for 
verification of safety properties where the goal is to prove that trajectories of a system starting 
from a given set of initial states will never reach a specified set of unsafe states [Alullj . This can 
be done by exhibiting an invariant set, containing the set of initial states, and whose intersection 
with the set of unsafe states is empty. For polynomial dynamics, these invariants are often given 
by semi-algebraic sets jPJP071 IPC081 ISanlOj . However, when the computation of an invariant set is 
part of a bigger process such as controller synthesis or safety verification, it is sometimes preferable 
to have invariants given by polytopes that are easier to manipulate [ALBHOTI ISDlOS] . For instance, 
for specific classes of polynomial systems such as multi-affine or quasi multi-affine systems, methods 
to obtain rectangular invariants have been developed in ^BH06l IATS09j . 

In this paper, we deal with the computation of polytopic invariant sets of polynomial dynamical 
systems. Let us remark that rectangles form a subclass of polytopic invariants and in some sense, 
our work extends the work of |BH06| IATS09j . More precisely, we shall consider a dynamical system 



1. Introduction 



of the form: 




x{t) = f{x{t)), x{t) E R' 



n 
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where / : M" — >• M" is a polynomial vector field. We consider the dynamics of ( 1.1 ) only on a bounded 
rectangle R of the state space MP; given a bounded polytope P <^ R with a set of facets {Fk\ k £ K}, 
it follows from the standard characterization of invariant sets (see e.g. |Aub91j ) that P is invariant 



for the dynamical system (1.1) if and only if 

(1.2) Wk £ K,yx£ Ffc, Qk ■ f{x) < 

where is the normal vector to Ff^. pointing outside P. As pointed out in |ATS09j and by application 
of Tarski's Theorem |Tar48j . this a decidable problem. However, the complexity of the decision 



procedure gives little hope for practical application. Let us remark that (1.2) can be reformulated 
as follows: 

(1.3) \/k G K, min -a^ • f{x) > 0. 

This consists in showing that the minimal values of the multivariate polynomials —at-f on the 
bounded polytopes are positive. Hence, if we are able to compute non-negative certified lower 
bounds of these minimal values, it is sufficient to prove that the polytope P is invariant for the 



dynamical system (1.1). 



In this paper, we establish linear programming (LP) relaxations of the optimization problems in 



(1.3). The main tool we use is the blossoming principle (see e.g. |Sei93j ) that essentially maps the 
set of polynomials to the set of symmetric multi-affine functions. The blossoming principle together 
with properties of multi-affine maps on bounded rectangles allows us to derive linear programs 



which makes it possible to compute lower bounds of the minimal values in (1.3). Our approach 
is conservative (the lower bound is not tight) but it is effective and may be sufficient for proving 
invariance of a polytope. Additionally, we will show how one can iteratively compute a polytopic 
invariant set using sensitivity analysis of linear programs. Finally, we show using a set of examples 
borrowed from biological applications, that our approach is effective in practice. 



2. Preliminaries 

In this section, we introduce notations and preliminary results that will be useful for subsequent 
discussions. 

2.1. Multi-afBne functions. Multi-affine functions form a particular class of multivariate polyno- 
mials. Essentially, a multi-affine function is a function which is affine in each of its variables when 
the other variables are regarded as constant: 

Definition 2.1. A multi-affine function p : M" — )• M is a multivariate polynomial in the variables 
xi, . . . ,Xn where the degree of p in each of the variable is at most 1. For x = {xi, . . . , Xn), 

P{x)= 

(«i,...,/„)e{o,i}" 
where pi,,„j„ e M for all (/i, ...,/„) G {0, 1}". 

Let R = Y[k=ii!!!^kJ ^k] be a rectangle of M", with X). < x^, for all A; G {1, . . . , n}; the set of vertices of 
R is V = Y[k=i{^kJ^k}- It is shown in |BH06j that a multi-affine function p is uniquely determined 
by its values at the vertices of a rectangle R. Moreover, for all x € R, p{x) is a convex combination 
of the values at the vertices, that is p{R) C CH{{p{v)\ v S V}) where CH{S) denotes the convex 
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hull of the set S. Let us remark that generally p{R) is not convex and P{R) / CH{{p{v)\ v G V}). 
Then, we have the following result: 

Lemma 2.2. Let p be a multi-affine function and R a rectangle with set of vertices V , then 
minp(x) = m.m.p{v). 

The previous result shows that optimizating a multi-affine function p over a rectangle R only requires 
finding the minimal value of p at the vertices of R. 

2.2. Blossoming principle. Let p : M" — )■ M be an arbitrary multivariate polynomial function, let 
6i, . . . ,5n denote the degree of p in the variables xi, . . . ,Xn respectively. Let A = {0,...,(5i}x---x 
{0, . . . , 5n}, then p{x) can be written under the form: 

P{x)= Ph,-,i„Xi ■ ■ ■ Xn 

(h,...,i„)eA 

where Pii^...^i„ G K for all (/i, . . . , /„) E A. We now present the blossoming principle which consists 
in mapping polynomials to symmetric multi-affine functions. Blossoms have been developed in the 
area of computer aided geometric design where they have numerous applications, most notably for 
spline curves and surfaces. As for our problem, the blossoming principle allows us to recast the 
optimization of polynomial functions as the optimization of a multi-affine functions for which we 



can use the fundamental property presented in Lemma 2.2 All the results in this section are quite 



standard (see [Sei93 j and references therein) and are therefore stated without proofs. 

Definition 2.3. The blossom or polar form of the polynomial p : M" — t- M is the function q : 

i=n 

Qiz)= Ph,...,inYlBi^,sM,i^---^Zi^sJ 

(«i,...,«„)eA i=l 

with 

Bl^sizi, . . . , Zs) = jyT Y 

where C{1,6) denotes the set of combinations of / elements in {1, ... ,6}. 

An example may help to understand the definition; the blossom of the polynomial p{x) = 3xi + 
2x2 + ^1^2 is 

q{z) = §(^1,1 + 21,2) + 2z2,i2;2,2^;2,3 

+ 5^1,l2l,2(22,l2;2,2 + 2:2,1^2,3 + Z2,2Z2,3)- 

We define a relation on M'^i+'''+'^": for z, z' G ]R'^i+'''+'^", with z = {zi^i, . . . , zi^s-^, ■■■ , Zn,i, ■■■ , Zn,5„) 
and z' = {z'l^i, . . . , z[ . . . , z^^^, . . ■ , z'^ ^^J, we denote z = z' if, for all k = 1, . . . , n, there exists 
a permutation tt^ such that {z^^i, . . . , z^^s^) = '^k{z'k i, ■ ■ ■ , z'f^ ^^). It is easy to see that 



equivalence relation. A characterization of blossoms that is equivalent to Definition 2.3 is given by 
the following proposition: 

Proposition 2.4. q : M^^^ — )• M is the blossom of the polynomial p : M" — )• M z/ and only if: 
(1) q is a multi-affine function; 
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(2) q is a symmetric function of its arguments: 

Vz ^ z', q{z) = q{z'); 

(3) q satisfies the diagonal property: 

q{zi, . . . , Zi, . . . , Zn, . . . , Zn) = p{zi, . . . , Zn). 

Let R = Yl^=ii^k'^k] be a rectangle of M", with j;^ < x^, for all k S {!,..., n}; we define 
the associated rectangle of M^^"' defined as R' = Yit—ii^k^^k]^'' ^^'^ its set of vertices V = 
Y\k=i{xi:,XkY''. For V = (-ui,!, . . . . . .,Vn,i, ■ ■ ■,Vn,sJ ^ V and k G {1, . . . ,n}, lk{v) denotes 

the number of elements Vk,i, ■ ■ ■ ,Vk,Sk that are equal to Xk- It is easy to verify that for v,v' G V , 

V ^ v' if and only if k{v) = lk{v') for all A; G {1, . . . , n}. We denote by v' = {V / ^) the set of 
equivalence classes of the relation = on the set V] V has (5i + 1) x • • • x (5„ + 1) elements. 

From the previous discussion, for k G {!,..., n} and v G V' , lk{v) has the same value for all 

V € V, with a sl ight abuse of notation we denote this value hiv)- Also from the second property in 
Proposition |2.4| for v G v' , q{v) has the same value for all v £v, let us denote this value q{v). 

Proposition 2.5. The values q{v) forv G V' are the coordinates of the polynomial p in the Bernstein 
basis: 

k=n 

veV' ^=1 
where Bi^s ^'"e the Bernstein polynomials 

Bi^s{y) = {\)y\i-vf~\ Zg{i,...,<5}. 

The previous result can be useful when one needs to compute the values qiv) for v G V' . The 
explicit computation of the blossom q{z) (which can count up to 2^^"' terms) is computationally 
expensive. However, Proposition |2.5| states that it is sufficient to compute the ((^i + 1) x • • • x [6n + 1) 
coordinates of the polynomial p{x) in the Bernstein basis. This can be done simply by solving a 
system of linear equations. 




3. LP Relaxations for Optimization of Polynomials on Polytopes 



As stated in the introduction, the verification of polytopic invariants for polynomial dynamical 
systems can be handled by solving a set of problems of optimization of multivariate polynomials on 
bounded polytopes given by (1.3). Therefore, in this section, we consider the following problem: 

min p{x) 

, over X G -R, 

under ai ■ x < bi, i ^ I, 
Cj ■ X = dj , j £ J- 

where p : M" — )• M is a multivariate polynomial, i? is a rectangle of with set of vertices V; 
/ = {!,..., mj} and J = {1, . . . , mj} are sets of indices; Oj G M", bi G M, for all z G I and Cj G M", 
dj G M, for all j G J. Let us remark that even though the polytope defined by the constraints indexed 
by I and J is unbounded in M", the fact that we consider x £ R which is a bounded rectangle of 
results in an optimization problem on a bounded (not necessarily full dimensional) polytope of 
M". We will assume that the problem is feasible: there exists x £ R satisfying all the constraints. 
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Problems defined in (1.3) can be recasted under the form (3.1) with polynomial p{x) = — • /(x) 
and linear inequality and equality constraints indexed by / and J describing the facet -F^. 

As the function p is usually non-convex, this may be a non-trivial problem to solve. Let us remark 
that as far as verification of polytopic invariants is concerned, we are not interested in computing 



the solution of problem (3.1) (i.e. x* ^ R satisfying constraints and minimizing p). Indeed, it is 
sufficient to compute the optimal value of (3.1), that is p* = p{x*), or at least a certified lower bound 



of the optimal value. We will first show how this can be done if p belongs to the particular class of 
multi-affine functions. Then, we will extend these results to arbitrary polynomial functions. 



3.1. Optimization of multi-afHne functions. In the following, we show how to compute, using 



linear programming, a certified lower bounded of the optimal value p* of (3.1) where p is a multi- 



affine function. The linear program is derived through Lagrangian duality. We start by writing the 

fJ,j{Cj ■ X 



Lagrangian of problem (3.1): 



A, //) = p{x) + ^ \i{ai ■ X - bi) + 

for all j G J. Then, the dual formulation of 



i6J 



where x G i?, Aj > for all i ^ I, and /ij G 
problem (3.1 ) is 



max 



(3.2) 



minL(x, A, /x) 

over A G W^' , fi G 
under Aj > 0, 



pmj 



i G /. 



Since (3.1) is feasible, the optimal value of (3.2) is bounded, it is denoted d* G M. It is well-known 



from duality theory (see e.g. [BV04j) that we have d* < p* . A multi-affine function is generally 
non-convex; then, we cannot expect strong duality (i.e. d* = p*) in general. The following result 
shows that problem (|3.2|) can be recasted as a linear program: 



Proposition 3.1. Let p 

to the linear program: 



(3.3) 



he a multi-affine function, the dual of problem (3.1) is equivalent 



, n G 



max t 

over t G M, AG MT' 
under Xi > 0, 

t < p{v) + ^ Aj(a 



pmj 



bi) 



dj), 



i G I, 



v£V. 



Proof. Let us remark that the Lagrangian L{x, A, /u) is a multi-affine function of x. Then, it follows 



from Lemma 2.2 that the minimum of L(x, X, /j.) over x G i? is reached at one of vertices of R: 



minL(x, A, fi) = m.mL{v, A, n). 

x&R v^V 



Then, problem (3.2) consist in optimizing a piecewise linear function under a set of linear constraints 
which it is straightforward (see [BV04j . pp 150-151) to formulate as the linear program (3.3). □ 



We have presented a simple approach to compute a certified lower bound of the solution p* of (3.1 ). 
Though conservative, our approach relies on linear programming and is therefore effective. 
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3.2. Optimization of polynomial functions. We consider problem (3.1 ) where p : 



IS now 



an arbitrary multivariate polynomial function. We use the blossoming principle to define a problem 



equivalent to (3.1) and involving q the blossom of p. We use the same notations as in section 2.2 



Then, from the third property in Proposition 2.4 problem (3.1) is equivalent to 



(3.4) 



mm 
over 
under 



q[z) 
z G R', 
CLi ■ z < bi 
c/ ■ z = d 



e-ki ■ z 



'J' 
0, 



j G J. 

A; G {1, . . . 

ZE{l,...,5fc-l}. 



where a'^ = 

and ek,i G 
/ G {1,.. 



^ Si ■ 



<5i ' 



(SiA 
V <5i • 



Si 



-Zkj+i,ior allz G 



Sn ' 



,1f)>foriG J; 
%/!:G{l,...,n}, 



P<5iH h<5n g^j,g ^jjg vectors such that ek^vz = zj^^i 
.Sk-l}. 

Now, since the blossom of a multivariate polynomial is a multi-affine function, we can remark that 
problem (3.4) is similar to those considered in Section 3.1 Then, we can use Proposition 3.1 to 



obtain its dual, which is given by the following linear program: 

max t 

over t£R, XgW^', fi£W^J, 
«gM('5i-l)+-+{5n-l), 

under A, > 0, 

t<q{v) + J2^i(.4-v-bi) 

fce{l,...,n} i6{l,...,5fc-l} 



i G /, 



(3.5) 



V G v. 



By Proposition 3.1 the optimal value of this linear program gives a certified lower bound of the 



optimal value of the polynomial optimization problem (3.1). However, we shall not solve directly 



the linear program (3.5) as the equivalence relation = defined in section 2.2 can be used to exploit 



the symmetries and reduce the complexity of the linear program. Let us remark that for all v 



i £ I and j G J, a'- 



a' • v' and c'- • v 



■ V . 



Then, for v € V , a'^- v (respectively c'- • v) has the 



same value for all v G v, let us denote this value a[ ■ v (respectively c'- ■ v). 



Theorem 3.2. Let p : 



be a polynomial and q : 



»(5lH \-Sn 



its blossom. The optimal 



value of the linear program (3.5) is equal to the optimal value d* of: 



(3.6) 



max t 

over t£R, XgW^', fiGW^J, 
under Aj > 0, i G I, 

t < q{v) + ^ Ai(a- -v-bi) 



Moreover, d* < p* where p* is the optimal value of (3.1). 
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of constraints indexed by G in (3.5 ) is 2^^ 



The proof is stated in appendix. Let us highlight the gain of solving (3.6) in place of (3.5) 

the decision variables a^i in problem (3.5) do not appear anymore in (3.6) 

■+<5, 



Firstly, 

Secondly, the number 
whereas the number of constraints indexed by 
V ^V' \s only (5i + 1) x • • • x ((5„ + 1). In fact, the linear program (3.5) has m/ + mj + {5i — 1) + 
■ ■ ■ + {^n — 1) + 1 variables and 2*^^^ + mi inequality constraints whereas the linear program (3.6 ) 
has only m/ + mj + 1 variables and ((5i + 1) x • • • x [5n + 1) + mi inequality constraints. As for 
the computation of the values in order to keep the computational cost as lows as possible, 

one should avoid computing explicitly the blossom of p; it is better to use the method suggested by 



Proposition 2.5 This way, the overall cost of the optimization procedure will remain polynomial in 
the degrees of p, though exponential in the dimension of the state space M". 



3.3. Sensitivity analysis. An interesting feature of Lagrangian duality is that it enables sensitivity 
analysis (see e.g. |B V04j ) . In this section, we are interested in analyzing the variations of the optimal 
value of (3.1), or of its lower bound, under modifications of the polytope. This will be used in the 
next section for the computation of polytopic invariants for polynomial dynamical systems. More 



precisely, we consider the following variation of problem (3.1): 

p[x) 



mm 
over 
under 



where G M, for all i G / and = 
problem (3.1) when a = and /3 



ai- X <hi + Ui, 
Cj ' X — dj I ^ 

(3j G M for aU j G J. 
: 0. We assume that problem (3.7) is feasible as well. Let p* and 



i G /, 

j e J, 

This problem coincides with the original 



P*(q,/3) denote the optimal values of problems (3.1) and (3.7), respectively. Let d* and d*{a,f3) be 
the lower bounds of p* and p*{a,(3) obtained by application of Theorem 3.2 The following result 



shows how the solution of (3.6) allows us to compute a lower bound of d*{a, f3) and thus of p*{a, (3). 



Theorem 3.3. Let d* and {t*,X*,fi*) be the optimal value and an optimal solution of the linear 
program (g.6| j. Then, for all a G M."^' and j3 G W^^-' , such that (3.1) is feasible we have: p*{a,j3) > 
d*{a,f3) >d* - X* - a- ■ 13. 



Proof. By applying Theorem 3.2 to (3.7), we have that is the optimal value of 



(3.8) 



max t 

over t G M, A G M™^ ;U G 

under Aj > 0, i £ I, 

t < q{v) + ^ Aj(a- - v-bi-Oi) 



The fact that p*{a,P) > d*{a,f3) is a consequence of Theorem 3.2 Let {t*,X*,ii*) be an optimal 



solution of the problem (|3.6[) , one can verify easily that (t * — A* • a — /U* • /? , A* , /x* ) is feasible for (3.8). 



It follows that d*{a, f3) >t* — X* ■ a — n* ■ f3 which leads to the expected inequality since d* = t* . □ 



3.4. Examples and comparison. Before using the method desribed previously for the verifica- 
tion and the computation of polytopic invariants for polynomial dynamical systems, we provide a 
brief comparison with two existing relaxation methods for polynomial optimization: the first one is 
the Reformulation Linearisation Technique (RLT) introduced by Sherali in |ST9H IST97] which is 
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also based on linear programming; the second one was introduced by Lasserre in [LasOlj and uses 
ralexations uner the form of Linear Matrix Inequalities (LMI) which are solved using semi definite 
programming. We compare the three methods using two different examples. 

We first consider the following constrained 3-dimensional problem |ST91] : 

min X1X2X3 + — 2xiX2 — Sxix^ 

+5X23^3 - + 5X2 + X3 

(3.9) over x = {xi,X2,X3) G [2,5] x [0,10] x [4,8], 

under 4xi + 8x2 + X3 < 20, 
xi + 2x2 + a^s > 1. 

Optimal value and solution of this problem are p* = —119 and x* = (xi*, X2*, X3*) = (3,0,8). 
Characteristics of the three relaxtions methods are collected in the table in Figure [T] 





RLT 


LMI 


Blossom 


Constraints 


56+2 


10x10+8(4x4) 


18+2 


Variables 


19 


34 


3 


Optimal value 


-120 


-119 


-120 


CPU time (sec) 


0.049 


0.548 


0.024 



Figure 1. Characteristics of the three methods for problem (3.9). 



We can see that our approach is the one leading to the simplest relaxation (problem with fewer 
variables and constraints). The time needed to compute a lower bound of the optimal value is 
therefore less than for the other approaches. The computed lower bound is the same than that 
computed by the RLT method but not as good as that computed by the LMI method which finds 
the exact minimum. Though, on that example, the gap between the computed lower bounds by the 
three methods is small. This is unfortunately not always the case as shown on the next example. 

Let us consider the following unconstrained 1-dimensional problem |VF92j (see also |ST97j ): 

min - 3y3 - 1.5^^ + lOy 
over y € [—5, 5], 

Optimal value and solution are p* = —7.5 and y* = —1. Characteristics of the three relaxation 
methods are collected in the table in Figure [2| 





RLT 


LMI 


Blossom 


Constraints 


5 


3x3+2(2x2) 


5 


Variables 


4 


4 


1 


Optimal value 


-908.3 


-7.5 


-837.5 


CPU time (sec) 


0.014 


0.246 


0.015 



Figure 2. Characterisitcs of the three methods for problem (3.10) 



The lower bounds provided by the two methods based on linear programming are far from the optimal 
value. Notice that for our approach this huge gap is typical from problems where the optimal solution 
is far from the boundary of the optimization domain. In that case, coupling our approach with a 
domain decomposition method based on branch and bound as suggested in |ST91] would strongly 
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improve our optimal value. In all cases, the LMI method gives better lower bounds but needs more 
computational resources. As for the problem of computing invariants for polynomial systems, it 
appears, as shown in Section [5| that our linear relaxations are often sufficient. 

4. Computation of Polytopic Invariants for Polynomial Systems 

In the following, we show how the results developed in the previous section can be used for the 
computation of polytopic invariants for polynomial systems. Let us consider the dynamical system 



(1.1), and let R = [xi,xi] x • • • x with < x^ for all k G {1, . . . , n} be a rectangle of M", 

delimiting a region of interest for studying the dynamics. Our goal is to compute a polytope P C i? 



invariant for (1.1). To restrict the search space, we shall use parametrized template expressions for 
P. We will impose the orientation of the facets of polytope P by choosing normal vectors in the set 
{ak G M*^! k € where K = {1,. . . ,mK} is a set of indices. Then, polytope P can be written 
under the form 

P = {x G M"| afc • X < \/k G K} 
where the vector b £ M™^, to be determined, specifies the position of the facets. The facets of P are 
denoted by for k £ K, where 

Fk = {x £ M^^l ttk • X = and ai ■ x < bi, yi £ K \ {k}}. 

The proposed approach for the computation of an invariant is iterative. Each iteration consists of two 
main steps. First, we try to verify that the polytope P is invariant for the dynamical system ( |1.1[ ). 
If we fail to verify that P is an invariant, we use sensitivity analysis to modify the vector b (and thus 
P) in order to find an invariant polytope. 

4.1. Polytopic invariant verification. As stated in the introduction, P is an invariant set of the 



dynamical system (1.1) if and only if 



yk £ K, min —ak ■ f{x) > 0. 

We assume that all the facets Fk are not empty. Since for all k £ K, Q P Q R, then this 
problem is equivalent to showing that the optimal values p^, of the following optimization problems 
are non- negative for all k £ K: 

min -ak ■ f{x) 
,^ over X £ R, 

under ai ■ x < bi, i £ K \ {k}, 
a-k - x = bk- 



Since —ak • / is a multivariate polynomial, this problem is similar to (3.1 ). Therefore, by application 
of Theorem |3.2[ we have the following result: 



Proposition 4.1. For k £ K , let qk be the blossom of the multivariate polynomial -Ok ■ f and let 

d*^ be the optimal value of the linear program: 

max t 

over t G M, AG Mr'< , 
(4.2) under Aj > 0, i£K\{k}, 

t<qk{v) + ^\i{ai-v-bi) v£V'. 



If for all k £ K, d*^>Q, then P is an invariant polytope for dynamical system (1.1). 
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Proof. By applying Theorem |3.2| to problem (4.1), we obtain that p^. > d^, for all k ^ K. Then 



> implies that p^, > 0, for all k ^ K, and therefore P is an invariant polytope. □ 
Remark 4.2. The degrees of the multivariat e po lynomials — • / may not be all the same. This 



results in vectors a^, and sets v' in problem (4.2) that depend on the index k £ K. It is possible to 



avoid this by defining as — • g where g is the blossom of the polynomial vector field / defined 
as follows. For i,j G {l,...,n}, let 5ij be the degree of Xj in the multivariate polynomial fj{x) 
and let 5i = maxjg|x_...^„} 5ij. It is possible to regard fj as a multivariate polynomial with degrees 
5i , . . . , (5„ poss ibly with some zero coefficients and to define the associated blossom gj as defined in 



Definition 2.3 Then, for j S {1, . . . , n}, gj are the components of the blossom g : M^^"' — M" of 
the polynomial vector field /. For k ^ K, = —aj. ■ g are multi-affine functions defined on M^^^ 



with similar properties to the blossom of — • / that can be used in problem (4.2). 



Invariance of a polytope P can be verified by solving a set of linear programs (one per facet of P). 
In the case when we fail to verify that the polytope is invariant, sensitivity analysis may help us in 



modifying it in order to find an invariant polytope for (1.1). 



4.2. Polytope modification using sensitivity analysis. The verification of the invariance of the 
polytope P fails if < 0, for some k £ K. In that case, we would like to know how to modify the 
vector b (and thus P) in order to find an invariant polytope. For a G i?™^, let Pa be the polytope 
given by 

Pa = {x e M"| Uk- X <bk + ak, V/c G K}. 



For a = 0, we recover the polytope P, we would like to find a such that Pa is an invariant for ( 1.1 ) 
We impose additional constraints on Pa'- 

• constraints of form 6fc + Ofc < bk ensures Pa ^ R; 

• constraints of form bk <bk + ak ensures Pa 7^ 0; 

• —£ < Ok < £ ensures that Pa is close to P, where e is a parameter that can be tuned. 



Denoting for k G K, dtia) the optimal values of problems (4.2) for the polytope Pa, the sensitivity 



analysis in Theorem 3.3 gives us (i|.(a) > dt + A|, • a, for all k £ K. where d*^ and (i^, A^) are the 



optimal values and solutions of problems (4.2) for polytope P and k ^ K. Then, by Proposition 4.1 



for Pa to be an invariant polytope for dynamical system (1.1), it is sufficient that for all k ^ K 



+ A^ • a > 0. In order to find a suitable a, we can solve the following problem: 

max min (dl + At • a) 

over aGM^-ff, 

under < < Ofc, k £ K 

where = max(— e, bk — bk) and Ok = min(e, bk — bk)- This problem can be recasted as the following 
linear program: 

max t 

, . over t G M, a G M*^^', 

^ ' under t < d*^ - Xk* ■ a, k £ K, 

Ok < Oik < Ok, k £ K. 

Let {t*,a*) be an optimal solution of this linear program. If the optimal value is non-negative then 



it is sufficient to prove that Pa* is an invariant for the dynamical system (1.1). If the optimal value 



is strictly negative, then we go back to the verification stage with P = Pa* and iterate the approach. 
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Remark 4.3. Let us remark that the polytope Pa* computed by solving (4.3) may have empty facets. 
This results, for the empty facet F^, in an unbounded value d1{a*) = +oo. In order to avoid such 
situations, it is useful to replace a* by a* such that P^* has no empty facet and P^* = Pa* (see 
Figure [3]). Again, this can be done by solving a set of linear programs. 




Figure 3. The polytope Pa* may have empty facets (center polytope), we replace 
a* by a* such that Pa* has no empty facet and Pa* = Pa* (right polytope). 



4.3. Related work. Computation of invariants for polynomial dynamical systems is often ap- 
proached using semi-definite programming via sum of squares relaxations. Most of the literature 
on the subject deals with the computation of Lyapunov functions, whose level sets are invariant (see 
e.g. I JW03t iiWLWOS] ) . A similar approach for the computation of invariant sets that do not contain 
any stable equilibrium point can be found in |PJP07j . In these works, the invariant sets are semi- 
algebraic sets described by polynomial inequalities. However, polytopes are easier to manipulate and 
as explained in [ALBHOT] , it is sometimes preferable to have invariants described by polytopes rather 
than semi-algebraic sets. In particular, polytopes with fixed facet directions given by a template has 
been shown very useful for the computation of invariant sets for hybrid systems with affine dynamics 
in |SDI0 8] and for multiaffine systems in jBH06|, lATSOOj . In some sense, our work builds on and 
extend these approaches to the class of polynomial dynamical systems. 

It is worth mentioning the work on polyhedral Lyapunov functions for the class of dynamical systems 
described by linear differential inclusions. For this class of systems, it can be shown that existence 
of an invariant set containing the origin is equivalent to Lyapunov stability |Bla99) . Methods relying 
on linear programming, for computing polyhedral invariants for linear differential inclusions have 
been developed (see e.g. |Bla95l [PolOOj). We would like to point out that our approach can be 
easily adapted for the computation of invariant sets for polynomial (and thus linear) differential 
inclusions. Hence, our approach can be seen as a generalization of the work mentioned above though 
the algorithms used for the computation of the polyhedral invariants are quite different. 



5. Examples 

We implemented our approach in Matlab; in the following, we show for a set of examples borrowed 
from biological applications, that our approach is effective in practice. All the reported computations 
take a few seconds. 

5.1. FitzHugh-Nagumo neuron model. We applied our approach to the FitzHugh-Nagumo 
model [Fit61j , a polynomial dynamical system modelling the electrical activity of a neuron: 

±1 = Xi - xl/3 - X2 + I, 

±2 = 0.08(xi +0.7-0.8x2), 
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where model parameter / is taken equal to |. This system is known to have a limit cycle. Using our 
approach, we synthesized an invariant polytope containing the limit cycle. Working in the rectangle 
[—2.5,2.5] X [—1.5,3.5], we found an invariant polytope with 8 facets with uniformly distributed 
orientations (see Figure [4]). Starting from the set represented in dashed line, our approach needs 15 
iterations to find the invariant polytope depicted in plain line. We can check on the figure that it is 
effectively an invariant. Let us remark that this invariant polytope P together with the existence of 
an unstable equilibrium inside P provides by application of the Poincare Bendixon theorem a formal 
proof of the existence of a limit cycle inside the polytope P. 

At each iteration, we need to solve for invariant verification, tjik = 8 linear programs of the form 



(4.2) with rriK + 1 = 9 variables and ttik — 1 + ((^i + 1) x (62 + 1) = 15 inequality constraints. For 
invariant synthesis using sensitivity analysis, we need to solve at each iteration 1 linear program with 
rriK + 1 = 9 variables and 2mK = 16 inequality constraints. 




Figure 4. Polytopic invariant for the FitzHugh-Nagumo model (represented in plain 
line) obtained after 15 iterations starting from the dashed polytope. The computed 
invariant contains the limit cycle. 



5.2. Phytoplankton growth model. We now consider a model of Phytoplankton growth [BG02| : 

XI = l-xi-3^, 
±2 = (2x3 - 1)3^2, 
X3 = ^-2xi. 

This system has a stable equilibrium. Using our approach, we synthesized an invariant polytope 
containing the equilibrium. Working in the rectangle [0,3] x [—0.1,2] x [0,0.6], we were able to find 
an invariant polytope with mx = 18 facets a regular octagon (see Figure [S]). Starting from the 
polytope represented in left part of the figure, our approach needs 11 iterations to find the invariant 
polytope depicted in the right part of the figure. We can check on the figure that it is indeed an 
invariant. 

At each iteration, we need to solve for invariant verification, mx = 18 linear programs of the form 



( 4.2 ) with rriK + 1 = 19 variables and mx — 1 + (5i + 1) x (^2 + 1) x ((^3 + 1) = 29 inequality constraints. 
For invariant synthesis using sensitivity analysis, we need to solve at each iteration 1 linear program 
with rriK + 1 = 19 variables and 2mK = 36 inequality constraints. 
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Figure 5. Polytopic invariant for the Phytoplankton growth model (on the right) 
obtained after 11 iterations starting from the polytope on the left. 

6. Conclusion 

In this paper, we have presented an approach based on linear programming for the computation of 
polytopic invariants for polynomial dynamical systems. It uses the blossoming principle for polyno- 
mials, properties of multi-affine functions, Lagrangian duality and sensitivity analysis. Though our 
approach is conservative (we may fail to verify invariance of a polytope), we have shown on several 
examples that it can be useful in practical applications. In the future, we plan to use similar ideas 
for control synthesis for polynomial dynamical systems. 
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Appendix - Proof of Theorem 13.21 



Let d\ be the optimal value of (3.5) and the optimal value of (3.6), we want to show that d\ = d!^. 
We start by remarking that (3.5) and (3.6) are linear programs therefore their optimal values are 
equal to that of their dual problems. One can show verify the dual of (3.5) is 



(6.1) 



min 






over 






under 


Vv > 0, 






-^VvV < hi, 


i G I, 










c'j ■ ^VvV = dj, 


je J, 




veV 






ek,i ■ '^VvV = 0, 


k£{l 
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Similarly, the dual of (3.6) is 



mm 



over 



Z £ ig('5i+l)x-x(5„+l) 



(6.2) 



under Zv >0, v £ V , 

Y Zy = l, 

veV' 

a'i ■Y^zjjv< hi, i e I, 

c'j ■ 'Y z^v = dj, j e J. 

tJeF' 

p2*^lH ^^ri 



We first show that < Let y G M be a feasible point for problem (6.1) such that 

= Y^veV Vvliv)- For v G v' , let z^ = YlvevVy clear that > 0. Further, 

12^^= Yl Y_y^ = X] = 1' 

a- • ^ z^t; = ^ ^ y„(a'i • ?j) = ^ y„(a'i • v) < bi, 



and 



Therefore, z is feasible for problem ( |6.2[ ). Finally, since for all f G u, = q{v), it follows that 

Therefore, < dl. We now show that d^ < d^. Let z G m(^i+i)x-x('5"+i) be a feasible point for 
problem (6.2) such that d2 = J2y,=v' ^^^'?(^)- ^(^) denote the number of vertices v Gv, then for 
all V £v, let yy = Zy/n{y). It is clear > and 

Y.y- = Y. Y.y- = Y.^^ = ^- 

vdV y^y' vdv ^gy' 

We also have that 

«i • X] = YYy"^^'i'^^=YY 

= Y ^"sia'i-v) <hi, 

and similarly we can show • ^^gy/ VvV = dj. Further, 

E\ ^ \ ^ Zy 

v&V y^y' v&v ^ ' 

\ ^ Zy I \ ^ 



Zv I I _N 
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By remarking, that for all v G v' , e^^i ■ Yl^fZyV = 0, it follows that Ck^i ■ ^^gy/ UvV = 0. Therefore, y 
is feasible for problem ( |6.1[ ). Finally, 



This proves that < ^2 and d* = d^ = The fact that d* < is a consequence of Proposition 3.1 
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